/*
  tempgraph
  A command line utility for plotting temperature data
  Copyright (C) 2009-2013 Bob Mottram
  bob@sluggish.dyndns.org

  This program is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include <iostream>
#include <stdio.h>
#include <vector>
#include <string>
#include <algorithm>

#include "anyoption.h"
#include "tempdata.h"
#include "gridcell.h"
#include "ghcn.h"
#include "gnuplot.h"
#include "globalgrid.h"

#define VERSION 1.40

#ifndef _WIN32
    #define TEMPGRAPH_BINARY "/usr/bin/tempgraph"
    #define COPY_FILE "cp"
    #define DELETE_FILE "rm -f"
    #define MAKE_DIR "mkdir"
    #define PATH_SEPARATOR "/"
#else
    #define TEMPGRAPH_BINARY "C:\tempgraph"
    #define COPY_FILE "copy"
    #define DELETE_FILE "del"
    #define MAKE_DIR "mkdir"
    #define PATH_SEPARATOR "\\"
#endif

using namespace std;

string string_to_lower(string str)
{
    const char* ch = str.c_str();
    string s = "";
    for (int i = 0; i < (int)str.size(); i++) {
		s += tolower(ch[i]);
	}
	return s;
}

void parse_string(string str, vector<string> &v)
{
    const char* ch = str.c_str();
    string s = "";
    for (int i = 0; i < (int)str.size(); i++) {
    	if (!((ch[i] == '|') || (ch[i] == ','))) {
    		if (!((ch[i] == ' ') && (s == ""))) {
    			s += tolower(ch[i]);
    		}
    	}
    	else {
    		if (s != "") {
    		    v.push_back(s);
    		    s = "";
    		}
    	}
    }
    if (s != "") {
    	if (s != "") {
    	    v.push_back(s);
    	}
    }
}

void parse_string(string str, vector<float> &v)
{
    const char* ch = str.c_str();
    string s = "";
    bool negate = false;
    for (int i = 0; i < (int)str.size(); i++) {
    	if (!((ch[i] == '|') || (ch[i] == ','))) {
    		if (!((ch[i] == ' ') && (s == ""))) {
    			if ((toupper(ch[i]) == 'S') || (toupper(ch[i]) == 'E')) negate = true;
    			if (!((toupper(ch[i]) == 'S') || (toupper(ch[i]) == 'E') ||
    				  (toupper(ch[i]) == 'N') || (toupper(ch[i]) == 'W'))) {
    			    s += tolower(ch[i]);
    			}
    		}
    	}
    	else {
    		if (s != "") {
    			if (!negate)
    		        v.push_back(atof(s.c_str()));
    			else
    				v.push_back(-atof(s.c_str()));
    		    s = "";
    		    negate = false;
    		}
    	}
    }
    if (s != "") {
    	if (s != "") {
    		if (!negate)
    	        v.push_back(atof(s.c_str()));
    		else
    			v.push_back(-atof(s.c_str()));
    	}
    }
}

void parse_latitude(string str, float &min, float &max)
{
    const char* ch = str.c_str();
    string s = "";
    int ctr = 0;
    bool negate = false;
    for (int i = 0; i < (int)str.size(); i++) {
    	if (!((ch[i] == '|') || (ch[i] == ','))) {
    		if (!((ch[i] == ' ') && (s == ""))) {
    			if ((toupper(ch[i]) == 'S') || (toupper(ch[i]) == 'E')) negate = true;
    			if (!((toupper(ch[i]) == 'S') || (toupper(ch[i]) == 'E') ||
    				  (toupper(ch[i]) == 'N') || (toupper(ch[i]) == 'W'))) {
    			    s += tolower(ch[i]);
    			}
    		}
    	}
    	else {
    		if (s != "") {
    			float v = 0;
    			if (!negate)
    		        v = atof(s.c_str());
    			else
    				v = -atof(s.c_str());
    			if (ctr == 0)
    				min = v;
    			else
    				max = v;
    		    s = "";
    		    negate = false;
    		    ctr++;
    		}
    	}
    }
    if (s != "") {
		float v = 0;
		if (!negate)
			v = atof(s.c_str());
		else
			v = -atof(s.c_str());
		if (ctr == 0)
			min = v;
		else
			max = v;
    }
}

void parse_view(string str, float &longitude, float &latitude)
{
    const char* ch = str.c_str();
    string s = "";
    int ctr = 0;
	int north=-1;
    bool negate = false;
    for (int i = 0; i < (int)str.size(); i++) {
    	if (!((ch[i] == '|') || (ch[i] == ','))) {
    		if (!((ch[i] == ' ') && (s == ""))) {
    			if ((toupper(ch[i]) == 'S') || (toupper(ch[i]) == 'E')) {
					negate = true;
				}
    			if (!((toupper(ch[i]) == 'S') || (toupper(ch[i]) == 'E') ||
    				  (toupper(ch[i]) == 'N') || (toupper(ch[i]) == 'W'))) {
    			    s += tolower(ch[i]);
    			}
    			if ((toupper(ch[i]) == 'S') || (toupper(ch[i]) == 'N')) {
					north=1;
				}
    			if ((toupper(ch[i]) == 'E') || (toupper(ch[i]) == 'W')) {
					north=0;
				}

    		}
    	}
    	else {
    		if (s != "") {
    			float v = 0;
    			if (!negate)
    		        v = atof(s.c_str());
    			else
    				v = -atof(s.c_str());
				if (north==-1) {
					if (ctr == 0)
						longitude = v;
					else
						latitude = v;
				}
				else {
					if (north==1) {
						latitude = v;
					}
					else {
						longitude = v;
					}
				}
    		    s = "";
    		    negate = false;
    		    ctr++;
				north=-1;
    		}
    	}
    }
    if (s != "") {
		float v = 0;
		if (!negate)
			v = atof(s.c_str());
		else
			v = -atof(s.c_str());
		if (north==-1) {
			if (ctr == 0)
				longitude = v;
			else
				latitude = v;
		}
		else {
			if (north==1) {
				latitude = v;
			}
			else {
				longitude = v;
			}
		}
    }
}

void parse_range(string str, float &min, float &max)
{
    const char* ch = str.c_str();
    string s = "";
    int ctr = 0;
    for (int i = 0; i < (int)str.size(); i++) {
    	if (!((ch[i] == '|') || (ch[i] == ','))) {
    		if (!((ch[i] == ' ') && (s == ""))) {
				s += tolower(ch[i]);
    		}
    	}
    	else {
    		if (s != "") {
    			float v = 0;
   		        v = atof(s.c_str());
    			if (ctr == 0)
    				min = v;
    			else
    				max = v;
    		    s = "";
    		    ctr++;
    		}
    	}
    }
    if (s != "") {
		float v = 0;
		v = atof(s.c_str());
		if (ctr == 0)
			min = v;
		else
			max = v;
    }
}

string parse_population_class(string population_class)
{
	population_class = string_to_lower(population_class);
	if ((population_class=="r") ||
		(population_class=="rural")) {
		return "R";
	}
	if ((population_class=="u") ||
		(population_class=="urban")) {
		return "U";
	}
	if ((population_class=="s") ||
		(population_class=="suburban") || (population_class=="suburb") ||
		(population_class=="suburbs")) {
		return "S";
	}
	return population_class;
}

string append_population_class_to_title(string title, string population_class)
{
	if (population_class=="R") {
		return "Rural " + title;
	}
	if (population_class=="S") {
		return "Suburban " + title;
	}
	if (population_class=="U") {
		return "Urban " + title;
	}
	return title;
}

string append_area_to_title(string title, vector<float> &area)
{
	if (area.size()==4) {
		if ((area[1] == 0.0f) && (area[3] == 0.0f)) {
			title += " between latitudes ";
			char str[10];
			sprintf(str, "%.1f", fabs(area[0]));
			if (area[0] >= 0) {
				title += string(str) + "N - ";
			}
			else {
				title += string(str) + "S - ";
			}
			sprintf(str, "%.1f", fabs(area[2]));
			if (area[2] >= 0) {
				title += string(str) + "N,";
			}
			else {
				title += string(str) + "S,";
			}
		}
		else {
			title += " within geographical area ";
			char str[10];
			sprintf(str, "%.1f", fabs(area[0]));
			if (area[0] >= 0) {
				title += string(str) + "N,";
			}
			else {
				title += string(str) + "S,";
			}
			sprintf(str, "%.1f", fabs(area[1]));
			if (area[1] >= 0) {
				title += string(str) + "W - ";
			}
			else {
				title += string(str) + "E - ";
			}

			sprintf(str, "%.1f", fabs(area[2]));
			if (area[2] >= 0) {
				title += string(str) + "N,";
			}
			else {
				title += string(str) + "S,";
			}
			sprintf(str, "%.1f", fabs(area[3]));
			if (area[3] >= 0) {
				title += string(str) + "W";
			}
			else {
				title += string(str) + "E";
			}
		}
	}

	return title;
}

string append_countries_to_title(string title,
								 vector<long long int> &country_codes,
								 vector<string> &countries)
{
	if (country_codes.size()>0) {
		title += " for ";
		for (int i = 0; i < (int)country_codes.size(); i++) {
			title += countries[i];
			if (i < (int)country_codes.size()-1) {
				if (i < (int)country_codes.size()-2) {
					title += ", ";
				}
				else {
					title += " and ";
				}
			}
		}
	}
	return title;
}

// TODO: The main function is quite hacky and could do with some refactoring

int main(int argc, char* argv[]) {
	vector<long long int> country_code;
	vector<string> country_name;

	AnyOption *opt = new AnyOption();

	// help
	opt->addUsage( "Example: " );
	opt->addUsage( "  tempgraph --start 1900 --end 2000 -- ghcnversion 2 --ghcn v2.mean --countrycodes v2.country.codes --countries uk --minmax" );
	opt->addUsage( " " );
	opt->addUsage( "Usage: " );
	opt->addUsage( "" );
	opt->addUsage( "     --start                Starting year for the plot");
	opt->addUsage( "     --end                  Ending year for the plot");
	opt->addUsage( "     --ghcn                 Load GHCN file");
	opt->addUsage( "     --wmo                  Load weather stations file or metadata file");
	opt->addUsage( "     --indent               Subtitle indent in the range 0.0-1.0");
	opt->addUsage( "     --vpos                 Subtitle vertical position in the range 0.0-1.0");
	opt->addUsage( "     --countries            List of countries to be plotted");
	opt->addUsage( "     --countrycodes         File containing country codes");
	opt->addUsage( "     --area                 Geographical area within which to plot temperatures");
	opt->addUsage( "     --stations             List of weather station numbers to be plotted");
	opt->addUsage( "     --minmax               Show min and max temperatures on the graph");
	opt->addUsage( "     --fitline              Show the best fit line");
	opt->addUsage( "     --width                Width of the image in pixels");
	opt->addUsage( "     --height               Height of the image in pixels");
	opt->addUsage( "     --latitudes            Plot within a min and max latitude band");
	opt->addUsage( "     --latitudetemps        Plot average temperatures at every latitude");
	opt->addUsage( "     --altitudetemps        Plot average temperatures at every altitude");
	opt->addUsage( "     --altitudes            Plot an altitude histogram");
	opt->addUsage( "     --avaltitude           Plot average altitude over time");
	opt->addUsage( "     --activestations       Plot number of active stations per year");
	opt->addUsage( "     --equatordist          Plot average station distance from equator per year");
	opt->addUsage( "     --distribution         Plot the distribution of temperatures");
	opt->addUsage( "     --interval             Time interval for distribution plots in years");
	opt->addUsage( "     --temprange            Show results within a range of temperature values");
	opt->addUsage( "     --monthrange           Show results within a range of months");
	opt->addUsage( "     --refyears             Define reference years for anomaly calculation");
	opt->addUsage( "     --anomalies            Plot anomalies");
	opt->addUsage( "     --seasonal             Plot seasonal anomalies");
	opt->addUsage( "     --runningaverage       Show running average anomalies");
	opt->addUsage( "     --map                  Plot a world map showing anomalies");
	opt->addUsage( "     --map3d                Plot a 3D world map showing anomalies");
	opt->addUsage( "     --showcentres          Show the centres of grid cells");
	opt->addUsage( "     --view                 Specify viewing longitude and latitude for 3D world map");
	opt->addUsage( "     --year                 Year used to plot world map");
	opt->addUsage( "     --popclass             Weather station population class");
	opt->addUsage( "     --save                 Filename to save the graph image as");
	opt->addUsage( "     --kml                  Save weather stations as KML for use with Google Earth");
	opt->addUsage( "     --ghcnversion          GHCN monthly data version, 2 or 3");
	opt->addUsage( "  -V --version              Show version number");
	opt->addUsage( "     --help                 Show help");
	opt->addUsage( "" );

	opt->setOption(  "start" );
	opt->setOption(  "end" );
	opt->setOption(  "year" );
	opt->setOption(  "ghcn" );
	opt->setOption(  "wmo" );
	opt->setOption(  "countries" );
	opt->setOption(  "countrycodes" );
	opt->setOption(  "stations" );
	opt->setOption(  "area" );
	opt->setOption(  "save" );
	opt->setOption(  "width" );
	opt->setOption(  "height" );
	opt->setOption(  "kml" );
	opt->setOption(  "activestations" );
	opt->setOption(  "equatordist" );
	opt->setOption(  "distribution" );
	opt->setOption(  "interval" );
	opt->setOption(  "latitudes" );
	opt->setOption(  "temprange" );
	opt->setOption(  "refyears" );
	opt->setOption(  "anomalies" );
	opt->setOption(  "seasonal" );
	opt->setOption(  "monthrange" );
	opt->setOption(  "map" );
	opt->setOption(  "view" );
	opt->setOption(  "map3d" );
	opt->setOption(  "ghcnversion" );
	opt->setOption(  "popclass" );
	opt->setOption(  "latitudetemps" );
	opt->setOption(  "altitudetemps" );
	opt->setOption(  "altitudes" );
	opt->setOption(  "avaltitude" );
	opt->setOption(  "indent" );
	opt->setOption(  "vpos" );
	opt->setFlag(  "version", 'V' );
	opt->setFlag(  "minmax" );
	opt->setFlag(  "fitline" );
	opt->setFlag(  "runningaverage" );
	opt->setFlag(  "showcentres" );
	opt->setFlag(  "help" );

	opt->processCommandArgs(argc, argv);

	if(( ! opt->hasOptions()) || opt->getFlag( "help" ))
		{
			// print usage if no options
			opt->printUsage();
			delete opt;
			return(0);
		}

	if( opt->getFlag( "version" ) || opt->getFlag( 'V' ) )
		{
			printf("Version %f\n", VERSION);
			delete opt;
			return(0);
		}

	bool best_fit_line = false;
	if( opt->getFlag( "fitline" ) ) {
		best_fit_line = true;
	}

	bool show_minmax = false;
	if( opt->getFlag( "minmax" ) ) {
		show_minmax = true;
	}

	bool show_cell_centres = false;
	if( opt->getFlag( "showcentres" ) ) {
		show_cell_centres = true;
	}

	bool show_running_average = false;
	if( opt->getFlag( "runningaverage" ) ) {
		show_running_average = true;
	}

	string average_altitude_filename = "";
	if( opt->getValue( "avaltitude" ) != NULL  ) {
		average_altitude_filename = opt->getValue("avaltitude");
		if (average_altitude_filename == "") average_altitude_filename = "average_altitudes.png";
	}

	string station_altitudes_filename = "";
	if( opt->getValue( "altitudes" ) != NULL  ) {
		station_altitudes_filename = opt->getValue("altitudes");
		if (station_altitudes_filename == "") station_altitudes_filename = "station_altitudes.png";
	}

	string latitude_temperatures_filename = "";
	if( opt->getValue( "latitudetemps" ) != NULL  ) {
		latitude_temperatures_filename = opt->getValue("latitudetemps");
		if (latitude_temperatures_filename == "") latitude_temperatures_filename = "latitude_temperatures.png";
	}

	string altitude_temperatures_filename = "";
	if( opt->getValue( "altitudetemps" ) != NULL  ) {
		altitude_temperatures_filename = opt->getValue("altitudetemps");
		if (altitude_temperatures_filename == "") altitude_temperatures_filename = "altitude_temperatures.png";
	}

	string active_stations_filename = "";
	if( opt->getValue( "activestations" ) != NULL  ) {
		active_stations_filename = opt->getValue("activestations");
		if (active_stations_filename == "") active_stations_filename = "graph_active_stations.png";
	}

	string equator_dist_filename = "";
	if( opt->getValue( "equatordist" ) != NULL  ) {
		equator_dist_filename = opt->getValue("equatordist");
		if (equator_dist_filename == "") equator_dist_filename = "graph_dist_to_equator.png";
	}

	string distribution_filename = "";
	if( opt->getValue( "distribution" ) != NULL  ) {
		distribution_filename = opt->getValue("distribution");
		if (distribution_filename == "") distribution_filename = "graph_distribution.png";
	}

	string ghcn_filename = "";
	float ghcn_version = 3.0f;
	string country_codes_filename = "";
	if( opt->getValue( "ghcn" ) != NULL  ) {
		ghcn_filename = opt->getValue("ghcn");

		if( opt->getValue( "countrycodes" ) != NULL  ) {
			country_codes_filename = opt->getValue("countrycodes");
			ghcn::load_country_codes(country_codes_filename, country_code, country_name);
		}
	}
	if( opt->getValue( "ghcnversion" ) != NULL  ) {
		ghcn_version = atof(opt->getValue("ghcnversion"));
	}

	// Subtitle indent
	float subtitle_indent = 0.34f;
	if( opt->getValue( "indent" ) != NULL  ) {
		subtitle_indent = atoi(opt->getValue("indent"));
	}

	// Subtitle vertical position
	float subtitle_vpos = 0.94f;
	if( opt->getValue( "vpos" ) != NULL  ) {
		subtitle_vpos = atoi(opt->getValue("vpos"));
	}

	int image_width = 1024;
	if( opt->getValue( "width" ) != NULL  ) {
		image_width = atoi(opt->getValue("width"));
	}
	int image_height = 768;
	if( opt->getValue( "height" ) != NULL  ) {
		image_height = atoi(opt->getValue("height"));
	}

	int distribution_interval = 0;
	if( opt->getValue( "interval" ) != NULL  ) {
		distribution_interval = atoi(opt->getValue("interval"));
	}

	vector<string> countries;
	if( opt->getValue( "countries" ) != NULL  ) {
		string countries_str = opt->getValue("countries");
		parse_string(countries_str, countries);
	}

	vector<float> area;
	if( opt->getValue( "area" ) != NULL  ) {
		string area_str = opt->getValue("area");
		parse_string(area_str, area);
		if (((int)area.size() > 0) && ((int)area.size() < 4)) {
			printf("Only %d out of 4 area coordinates were specified.  Did you miss out a comma?\n", (int)area.size());
			delete opt;
			return(0);
		}
	}

	if( opt->getValue( "latitudes" ) != NULL  ) {
		string latitude_str = opt->getValue("latitudes");
		float min_latitude=40, max_latitude=50;
		parse_latitude(latitude_str, min_latitude, max_latitude);
		area.clear();
		area.push_back(min_latitude);
		area.push_back(0.0f);
		area.push_back(max_latitude);
		area.push_back(0.0f);
	}

	string population_class = "";
	if( opt->getValue( "popclass" ) != NULL  ) {
		population_class = parse_population_class(opt->getValue("popclass"));
	}

	float min_temp=-100, max_temp=100;
	if( opt->getValue( "temprange" ) != NULL  ) {
		string temprange_str = opt->getValue("temprange");
		parse_range(temprange_str, min_temp, max_temp);
		printf("Temperature range: %.2f - %.2f\n", min_temp, max_temp);
	}

	float view_longitude=0, view_latitude=50;
	if( opt->getValue( "view" ) != NULL  ) {
		string view_str = opt->getValue("view");
		parse_view(view_str, view_longitude, view_latitude);
		printf("Viewing angle: %.1fW - %.1fN\n", view_longitude, view_latitude);
	}

	// Reference years used for anomaly calculation
	int reference_start_year=1961, reference_end_year=1990;
	if( opt->getValue( "refyears" ) != NULL  ) {
		string refyears_str = opt->getValue("refyears");
		float y1=1961,y2=1990;
		parse_range(refyears_str, y1,y2);
		reference_start_year = (int)y1;
		reference_end_year = (int)y2;
		printf("Reference period: %d - %d\n", reference_start_year, reference_end_year);
	}

	// month numbers in the range 1-12
	int start_month=0, end_month=12;
	const string monthstr[12] = {
		"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"
	};
	if( opt->getValue( "monthrange" ) != NULL  ) {
		string monthrange_str = opt->getValue("monthrange");
		float m1=0,m2=12;
		parse_range(monthrange_str, m1,m2);
		start_month = (int)m1-1;
		end_month = (int)m2-1;		
		printf("Month range: %s - %s\n", monthstr[start_month].c_str(), monthstr[end_month].c_str());
	}

	vector<string> stations;
	vector<long long int> station_numbers;
	if( opt->getValue( "stations" ) != NULL  ) {
		string stations_str = opt->getValue("stations");
		parse_string(stations_str, stations);
		for (int i = 0; i < (int)stations.size(); i++) {
			long long int number = (long long int)atof(stations[i].c_str());
			if ((ghcn_version>=3) &&
				(number<=99999)) {
				// If the 5 number station identifier is used then append three zeros
				number *= 1000;
			}
			station_numbers.push_back(number);
		}
	}

	string anomalies_filename = "";
	if( opt->getValue( "anomalies" ) != NULL  ) {
		anomalies_filename = opt->getValue("anomalies");
		if (anomalies_filename=="") anomalies_filename = "graph_anomalies.png";
	}

	string seasonal_anomalies_filename="";
	if( opt->getValue( "seasonal" ) != NULL  ) {
		seasonal_anomalies_filename = opt->getValue("seasonal");
		if (seasonal_anomalies_filename=="") seasonal_anomalies_filename = "graph_seasonal_anomalies.png";
		anomalies_filename = seasonal_anomalies_filename;
	}

	string kml_filename = "";
	if( opt->getValue( "kml" ) != NULL  ) {
		kml_filename = opt->getValue("kml");
	}

	int start_year = 1770;
	if( opt->getValue( "start" ) != NULL  ) {
		start_year = atoi(opt->getValue("start"));
	}

	int end_year = 2009;
	if( opt->getValue( "end" ) != NULL  ) {
		end_year = atoi(opt->getValue("end"));
	}

	int selected_year = 2009;
	if( opt->getValue( "year" ) != NULL  ) {
		selected_year = atoi(opt->getValue("year"));
		if (selected_year>end_year) end_year=selected_year;
		if (start_year>reference_start_year) start_year=reference_start_year;
		if (end_year<reference_end_year) end_year=reference_end_year;
	}

	string world_map_filename = "";
	bool threed = false;
	if( opt->getValue( "map" ) != NULL  ) {
		world_map_filename = opt->getValue("map");
		if (world_map_filename=="") world_map_filename = "map_anomalies.png";
	}
	if( opt->getValue( "map3d" ) != NULL  ) {
		world_map_filename = opt->getValue("map3d");
		if (world_map_filename=="") world_map_filename = "map_anomalies.png";
		threed=true;
	}

	vector<tempdata> data;
	vector<long long int> country_codes;

	vector<stationdata> weather_stations;
	std::string weather_stations_filename = "";
	if( opt->getValue( "wmo" ) != NULL  ) {
		weather_stations_filename = opt->getValue("wmo");

		if (kml_filename!="") {
			ghcn::get_country_codes(
									countries,
									country_code,
									country_name,
									country_codes);
		}

		ghcn::load(
				   ghcn_filename,
				   ghcn_version,
				   start_year,
				   end_year,
				   country_codes,
				   station_numbers,
				   data, 0);

		ghcn::load_stations(
							weather_stations_filename,
							ghcn_version,
							weather_stations,
							data,
							country_codes,
							min_temp,max_temp,
							start_month,end_month,
							population_class);

		if ((int)weather_stations.size() == 0) {
			printf("Warning: no weather stations could be found within %s\n", weather_stations_filename.c_str());
		}
		else {
			printf("%d weather stations loaded\n", (int)weather_stations.size());

			if (((int)area.size() == 4) &&
				(anomalies_filename=="") &&
				(world_map_filename=="") &&
				(distribution_filename=="") &&
				(active_stations_filename=="") &&
				(equator_dist_filename=="") &&
				(altitude_temperatures_filename=="") &&
				(latitude_temperatures_filename=="") &&
				(station_altitudes_filename=="") &&
				(average_altitude_filename=="")) {

				vector<stationdata> returned_stations;
				ghcn::get_stations_within_area(
											   area[0], area[1],
											   area[2], area[3],
											   weather_stations,
											   returned_stations);
				if ((int)returned_stations.size() == 0) {
					printf("No stations are located within the specified area\n");
				}
				else {
					printf("%d stations were located within the specified area\n", (int)returned_stations.size());

					stations.clear();
					station_numbers.clear();
					for (int i = 0; i < (int)returned_stations.size(); i++) {
						for (int j = 0; j < (int)weather_stations.size(); j++) {
							if (weather_stations[j].number == returned_stations[i].number) {
								station_numbers.push_back(weather_stations[j].number);
								stations.push_back(weather_stations[j].name);
								break;
							}
						}
					}

					if (kml_filename != "") {
						// plot only the stations inside the specified area
						vector<stationdata> new_weather_stations;
						for (int i = 0; i < (int)returned_stations.size(); i++) {
							for (int j = 0; j < (int)weather_stations.size(); j++) {
								if (weather_stations[j].number == returned_stations[i].number) {
									new_weather_stations.push_back(weather_stations[j]);
									break;
								}
							}
						}
						weather_stations = new_weather_stations;
					}
				}
			}

			if (kml_filename != "") {
				ghcn::save_stations(
									kml_filename,
									weather_stations);
				printf("Weather stations saved to %s\n", kml_filename.c_str());
			}
		}
	}

	string stations_image_filename = "graph_stations.png";
	string countries_image_filename = "graph_countries.png";
	string averages_image_filename = "graph_average.png";
	if( opt->getValue( "save" ) != NULL  ) {
		averages_image_filename = opt->getValue("save");
	}

	string stations_filename = "station_numbers.txt";
	string title = "";
	string subtitle = "";
	char str[20];
	if (kml_filename=="") {
		if (ghcn_filename != "") {
			ghcn::get_country_codes(
									countries,
									country_code,
									country_name,
									country_codes);

			if (((int)countries.size() > 0) && ((int)country_code.size() == 0)) {
				if (country_codes_filename == "") {
					printf("No country code filename was specified\nCheck the --countrycodes option\n");
				}
				else {
					printf("No country codes were loaded\nCheck that the file %s exists\n", country_codes_filename.c_str());
				}
			}
			else {
				printf("Loading temperature data...\n");

				if ((int)stations.size() == 0) {

					ghcn::load(
							   ghcn_filename,
							   ghcn_version,
							   start_year,
							   end_year,
							   country_codes,
							   station_numbers,
							   data, 0);

					if ((int)data.size() == 0) {
						printf("No temperature data was found matching the current criteria\n");
					}
					else {

						printf("%d temperature measurements loaded\n", (int)data.size());

						// plot by country

						subtitle = "source: GHCN http://www.ncdc.noaa.gov/ghcnm/";

						if (average_altitude_filename!="") {
							char titlestr[100];
							string title_average_altitudes = "";
							if (show_running_average) {
								title_average_altitudes = "Running ";
							}
							sprintf(titlestr,"Average Altitudes between %d and %d",start_year,end_year);
							title_average_altitudes += titlestr;

							title_average_altitudes = append_population_class_to_title(title_average_altitudes, population_class);

							if ((int)country_codes.size() > 0) {
								title_average_altitudes =
									append_countries_to_title(title_average_altitudes,
															  country_codes, countries);
							}

							title_average_altitudes = append_area_to_title(title_average_altitudes, area);

							gnuplot::plot_average_altitudes(
															"plot.gnu", "plot.dat",
															average_altitude_filename,
															title_average_altitudes,
															subtitle, subtitle_indent, subtitle_vpos,
															data,
															country_codes,
															weather_stations,
															start_year,
															end_year,
															min_temp, max_temp,
															population_class,
															area,
															show_running_average,
															image_width,
															image_height);
							printf("Average altitudes saved to %s\n", average_altitude_filename.c_str());							
						}
						else {

							if (station_altitudes_filename!="") {
								char titlestr[100];
								string title_station_altitudes = "";
								sprintf(titlestr,"Weather Station Altitudes Histogram between %d and %d",start_year,end_year);
								title_station_altitudes += titlestr;

								title_station_altitudes = append_population_class_to_title(title_station_altitudes, population_class);

								if ((int)country_codes.size() > 0) {
									title_station_altitudes =
										append_countries_to_title(title_station_altitudes,
																  country_codes, countries);
								}

								title_station_altitudes = append_area_to_title(title_station_altitudes, area);

								gnuplot::plot_station_altitudes(
																"plot.gnu", "plot.dat",
																station_altitudes_filename,
																title_station_altitudes,
																subtitle, subtitle_indent, subtitle_vpos,
																data,
																country_codes,
																weather_stations,
																start_year,
																end_year,
																min_temp, max_temp,
																population_class,
																area,
																image_width,
																image_height);
								printf("Station altitudes saved to %s\n", station_altitudes_filename.c_str());
							}
							else {
						
								if (altitude_temperatures_filename!="") {
									char titlestr[100];
									string title_altitude_temperatures="";
									sprintf(titlestr,"Average Altitude Temperatures between %d and %d",start_year,end_year);
									title_altitude_temperatures += titlestr;
									title_altitude_temperatures = append_population_class_to_title(title_altitude_temperatures, population_class);
								
									gnuplot::plot_altitude_temperature(
																	   "plot.gnu", "plot.dat",
																	   altitude_temperatures_filename,
																	   title_altitude_temperatures,
																	   subtitle, subtitle_indent, subtitle_vpos,
																	   data,weather_stations,
																	   start_year,end_year,
																	   population_class,
																	   show_minmax,
																	   image_width,
																	   image_height);
									printf("Altitude temperatures saved to %s\n", altitude_temperatures_filename.c_str());
								}
								else {


									if (latitude_temperatures_filename!="") {
										char titlestr[100];
										string title_latitude_temperatures="";
										sprintf(titlestr,"Average latitude temperatures between %d and %d",start_year,end_year);
										title_latitude_temperatures += titlestr;
										title_latitude_temperatures = append_population_class_to_title(title_latitude_temperatures, population_class);

										gnuplot::plot_latitude_temperature(
																		   "plot.gnu", "plot.dat",
																		   latitude_temperatures_filename,
																		   title_latitude_temperatures,
																		   subtitle, subtitle_indent, subtitle_vpos,
																		   data,weather_stations,
																		   start_year,end_year,
																		   population_class,
																		   show_minmax,
																		   image_width,
																		   image_height);
										printf("Latitude temperatures saved to %s\n", latitude_temperatures_filename.c_str());
									}
									else {

										if ((active_stations_filename != "") ||
											(equator_dist_filename != "")) {
											// Active stations plot *********************************************************
											string title_active_stations = "Number of active weather stations per year";
											string title_equator_dist = "Average distance of active weather stations to the equator";

											title_active_stations = append_population_class_to_title(title_active_stations, population_class);
											title_equator_dist = append_population_class_to_title(title_equator_dist, population_class);

											if ((int)country_codes.size() > 0) {
												title_active_stations =
													append_countries_to_title(title_active_stations,
																			  country_codes, countries);
												title_equator_dist =
													append_countries_to_title(title_equator_dist,
																			  country_codes, countries);
											}

											title_active_stations = append_area_to_title(title_active_stations, area);
											title_equator_dist = append_area_to_title(title_equator_dist, area);

											gnuplot::plot_active_stations(
																		  "plot.gnu", "plot.dat",
																		  active_stations_filename,
																		  equator_dist_filename,
																		  title_active_stations,
																		  title_equator_dist,
																		  subtitle, subtitle_indent, subtitle_vpos,
																		  data,
																		  country_codes,
																		  weather_stations,
																		  start_year,
																		  end_year,
																		  min_temp, max_temp,
																		  population_class,
																		  area,
																		  image_width,
																		  image_height);
											if (active_stations_filename != "") printf("Active stations saved to %s\n", active_stations_filename.c_str());
											if (equator_dist_filename != "") printf("Equator distances saved to %s\n", equator_dist_filename.c_str());
										}
										else {
											if (world_map_filename!="") {
												// World map ********************************************************************
												vector<gridcell> grid;
												char titlestr[80];
												// create the grid
												globalgrid::load(data,
																 weather_stations,
																 grid, area,
																 country_codes,
																 population_class);

												if (country_codes.size()==0) {
													sprintf(titlestr,"World Temperature Anomalies %d",selected_year);
												}
												else {
													sprintf(titlestr,"Temperature Anomalies %d",selected_year);
												}
												title="";
												title += titlestr;							

												title = append_population_class_to_title(title, population_class);

												title =
													append_countries_to_title(title, country_codes, countries);

												title = append_area_to_title(title, area);



												gnuplot::plot_world_map(
																		"plot.gnu", "plot.dat", "world.dat",
																		world_map_filename,
																		title, subtitle, subtitle_indent, subtitle_vpos,
																		grid,
																		selected_year,
																		reference_start_year,
																		reference_end_year,
																		threed,
																		view_longitude, view_latitude,
																		show_cell_centres,
																		image_width,
																		image_height);
											}
											else {
												if (anomalies_filename!="") {
													// Anomalies plot ***************************************************************
													vector<gridcell> grid;
													char titlestr[60];
													// create the grid
													globalgrid::load(data,
																	 weather_stations,
																	 grid, area,
																	 country_codes,
																	 population_class);

													if (show_running_average) {
														sprintf(titlestr,"Running Average Temperature Anomalies %d - %d",start_year,end_year);
													}
													else {
														sprintf(titlestr,"Temperature Anomalies %d - %d",start_year,end_year);
													}
													title="";
													title += titlestr;

													title = append_population_class_to_title(title, population_class);

													title =
														append_countries_to_title(title, country_codes, countries);

													title = append_area_to_title(title, area);

													if (seasonal_anomalies_filename=="") {
														gnuplot::plot_anomalies(
																				"plot.gnu", "plot.dat",
																				anomalies_filename,
																				title, subtitle, subtitle_indent, subtitle_vpos,
																				grid,
																				start_year, end_year,
																				reference_start_year,
																				reference_end_year,
																				show_minmax,
																				show_running_average,
																				best_fit_line,
																				image_width, image_height);
													}
													else {
														gnuplot::plot_seasonal_anomalies(
																						 "plot.gnu", "plot.dat",
																						 seasonal_anomalies_filename,
																						 title, subtitle,
																						 subtitle_indent, subtitle_vpos, grid,
																						 start_year, end_year,
																						 reference_start_year,
																						 reference_end_year,
																						 show_running_average,
																						 image_width, image_height);
													}
												}
												else {
													// Distribution plot ************************************************************
													if (distribution_filename!="") {
														string distribution_countries_filename="";
														if (countries.size() == 0) {
															if (area.size()==0) {
																title = "Global Temperature Distribution";
															}
															else {
																title = "Temperature Distribution";
															}
														}
														else {
															title = "Temperature Distributions";
														}
														if (end_year>start_year) {
															sprintf(str,"%d - %d",start_year,end_year);
															title += " " + std::string(str);
														}
														else {
															sprintf(str,"%d",start_year);
															title += " " + std::string(str);
														}

														title = append_population_class_to_title(title, population_class);

														title =
															append_countries_to_title(title, country_codes, countries);
									
														title = append_area_to_title(title, area);

														vector<gridcell> grid;
														// create the grid
														globalgrid::load(data,
																		 weather_stations,
																		 grid, area,
																		 country_codes,
																		 population_class);

														gnuplot::plot_distribution(
																				   "plot.gnu", "plot.dat",
																				   distribution_filename,
																				   title, subtitle,
																				   subtitle_indent, subtitle_vpos, grid,
																				   start_year, end_year,
																				   reference_start_year,
																				   reference_end_year,
																				   distribution_interval,
																				   image_width, image_height);

													}
													else {
														// Plot globally or by country *******************************************
														if (countries.size() == 0) {
															title = "Global Historical Temperatures";
														}
														else {
															if (countries.size() == 1) {
																title = "Historical Temperatures for " + countries[0];
															}
															else {
																title = "Historical Temperatures by Country";
															}
														}
														gnuplot::plot_by_country_annual(
																						"plot.gnu", "plot.dat",
																						countries_image_filename,
																						averages_image_filename,
																						stations_filename,
																						title, subtitle,
																						subtitle_indent, subtitle_vpos, data,
																						country_codes, countries,
																						start_year, end_year,
																						min_temp, max_temp,
																						show_minmax, best_fit_line,
																						image_width, image_height);

														printf("Country graphs saved to %s\n", countries_image_filename.c_str());
														printf("Average graphs saved to %s\n", averages_image_filename.c_str());
														printf("Station numbers saved to %s\n", stations_filename.c_str());
													}
												}
											}
										}
									}
								}
							}
						}
					}
				}
				else {
					// plot by station numbers ****************************************************

					ghcn::load(
							   ghcn_filename,
							   ghcn_version,
							   start_year,
							   end_year,
							   country_codes,
							   station_numbers,
							   data, 1);

					if ((int)data.size() == 0) {
						printf("No temperature data was found matching the current criteria\n");
					}
					else {

						printf("%d temperature measurements loaded\n", (int)data.size());

						if (stations.size() == 1) {
							title = "Historical Temperatures for Weather Station " + stations[0];
						}
						else {
							if ((int)area.size() == 4) {
								if ((area[1] == 0.0f) && (area[3] == 0.0f)) {
									stations_image_filename = "";
								}
								title = "Average Temperatures";
								title = append_area_to_title(title, area);
							}
							else {
								title = "Average Temperatures for Weather Stations ";
								for (int i = 0; i < (int)stations.size(); i++) {
									if (i < (int)stations.size()-1) {
										title += stations[i];
										if (i < (int)stations.size()-2) {
											title += ", ";
										}
										else {
											title += " and " + stations[(int)stations.size()-1];
										}
									}
								}
							}
						}

						subtitle = "source: GHCN http://www.ncdc.noaa.gov/ghcnm/";
						if (stations_image_filename != "") {
							gnuplot::plot_by_station_annual(
															"plot.gnu", "plot.dat",
															stations_image_filename,
															averages_image_filename,
															title, subtitle,
															subtitle_indent, subtitle_vpos, data,
															station_numbers, stations,
															start_year, end_year,
															show_minmax, best_fit_line,
															image_width, image_height);
						}
						else {
							// latitudes
							gnuplot::plot_annual_only(
													  "plot.gnu", "plot.dat",
													  averages_image_filename,
													  title, subtitle,
													  subtitle_indent, subtitle_vpos, data,
													  station_numbers, stations,
													  start_year, end_year,
													  show_minmax, best_fit_line,
													  image_width, image_height);
						}

						if (stations_image_filename != "") printf("Weather Station graphs saved to %s\n", stations_image_filename.c_str());
						printf("Average graphs saved to %s\n", averages_image_filename.c_str());
					}

				}
			}
		}
		else {
			printf("No GHCN file was specified\n");
		}
	}

	return 0;
}

